XeF2 1 - 300eV Multi-orb processing¶
PH
v2d 11/11/20
- Added simulated branching ratios of various types.
v2c 03/11/20
- Added expt. branching ratios comparison (summed).
v2b 27/10/20
- Updating plot styles to HV.
- Added branching ratios & more literature data.
- Added ion yield data.
v2 07/10/20 (updated and distributed 26/10/20)
- Running manually on Stimpy, ePSproc-v1.3-dev with multiJob class + RF updated experimental data analysis.
- (See N2O preliminary data analysis for template, http://localhost:8888/lab/tree/dev/ePSproc/classDev/ePSproc_multijob_class_tests_N2O_011020_Stimpy.ipynb)
v1 27/03/20
- Running manually on AntonJr, ePSprc-v1.2 env.
- Code adapted from recent multi-E chuncked job notebooks.
Setup¶
!hostname
!conda env list
import sys
import os
from pathlib import Path
import numpy as np
# import epsproc as ep
import xarray as xr
import matplotlib.pyplot as plt
from datetime import datetime as dt
timeString = dt.now()
# For module testing, include path to module here, otherwise use global installation
if sys.platform == "win32":
modPath = r'D:\code\github\ePSproc' # Win test machine
winFlag = True
else:
modPath = r'/home/femtolab/github/ePSproc/' # Linux test machine
winFlag = False
sys.path.append(modPath)
import epsproc as ep
# Plotters
from epsproc.plot import hvPlotters
# Multijob class dev code
from epsproc.classes.multiJob import ePSmultiJob
hvPlotters.setPlotters(width = 700)
# import bokeh
# import holoviews as hv
# hv.extension('bokeh')
Load data¶
# # Scan for subdirs, based on existing routine in getFiles()
fileBase = r'D:\projects\ePolyScat\xef2\XeF2_highRes_wf' # Data dir on Stimpy
data = ePSmultiJob(fileBase, verbose = 0)
keys = [4,5,6] # Set for 4d data only
data.scanFiles(keys = keys)
data.jobsSummary()
System properties¶
data.molSummary()
Orbitals 21 - 25 are the Xe(4d) core levels, with $A_{1g}$(SG), $E_{1g}$(PG) and $E_{2g}$(DG) components, where the classifcations are in full $D_{\infty h}$ and ePolyScat labels in parenthesis. The electronic structure was computed in Gamess (HF/SPK-QZP-rel), and converged to a bond-length of 1.9373A, and 4d energies as shown in the table above (~75 - 76eV). (See Molecular orbitals below for plots of the orbitals.)
Below: literature values. Comparing orb 21 to the lower 4d SO binding energy of 70.34 eV, the current computational results are off from the experimental binding energies by ~6.2 eV.
Source: Bancroft, G. M., P. A. Malmquist, S. Svensson, E. Basilier, U. Gelius, and K. Siegbahn. “Gas-Phase ESCA Studies of Valence and Core Levels in Xenon Difluoride and Xenon Tetrafluoride.” Inorganic Chemistry 17, no. 6 (June 1978): 1595–99. https://doi.org/10.1021/ic50184a040.
These values also compare well with the current experimental results - shown below from the file XeF2_4d_LV_115eV_binding_energy_extractions_with_atomic_Xe.png
Plot cross-sections and betas¶
These are from ePolyScat's getCro function, and are LF (unaligned ensemble) results. This provides a good, if general, overview. More detailed results - per orbital and symmetry - are given in the addendum. The structure/oscillations in some channels at higher energies remain to be explored/explained!
Overview (all symmetries, length gauge)¶
# Comparitive plot over datasets (all symmetries only)
Etype = 'Ehv' # Set for Eke or Ehv energy scale
Erange=[50, 300] # Plot range (full range if not passed to function below)
data.plotGetCroComp(pType='SIGMA', Etype = Etype, Erange = Erange)
# Comparative plot over datasets (all symmetries only)
data.plotGetCroComp(pType='BETA', Etype=Etype, Erange=Erange)
Branching ratios¶
# Stack XS data to new data structure
# NOTE: this is not quite correct, since it forces all data to same Ehv axis, but OK for a quick hack, and easy to pull out branching ratios
dsXS = xr.Dataset() # Set blank dataset, this is easier for stacking, probably
lText = []
for key in data.data.keys():
subset = data.data[key]['XS'].sel({'Sym':'All', 'XC':'SIGMA'}) # Set XS data, all syms only
# NOTE currently missing full dataset resolution for orb24, so try interp. (2eV not 1eV step size)
# Note dropna to ensure no NaNs, see http://xarray.pydata.org/en/latest/interpolation.html#interpolating-arrays-with-nan
if key == 'orb24':
subset = subset.dropna('Eke').interp(Eke = dsXS.Eke, method = 'cubic')
# dsXS[key] = subset.copy() # Set .copy() for safety here!
dsXS[data.data[key]['jobNotes']['orbLabel']] = subset.copy() # Set .copy() for safety here!
# Convert to Xarray & normalise
dsXS = dsXS.to_array().rename({'variable':'channel'}).squeeze()
dsXS['total'] = dsXS.sum('channel') # Sum over channels
# Normalise...
dsXS = dsXS/dsXS['total']
# Plot
Erange = [75, 250]
dsXS.swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L').plot.line(x='Ehv');
These are broadly similar over the different gauges. There's quite a bit of low-E structure, and things smooth out at higher-E, with ~2:1 ratio between SG:PG or DG channels.
A detailed comparison with expt. is made below, for either the SO summed case, or a full SO simulation - this currently looks reasonable, but is still a work in progress.
Experimental data¶
- Set 1: RF processed results, 02/10/20
- Set 2: DH ion yields, 27/10/20
# Load experimental data
dataPathExpt = Path(r'D:\projects\XeF2_Soleil_2019\RF_data_analysis_021020')
# Set data attribs in dict similar to ePS results structure
dataFiles = {}
dataFiles['SIGMA'] = {'Data':'XS', 'File':r'XeF2_Xe4d_4d32_4d52_cross_sec_all_photon_energies_02102020.dat',
'States':['$4d_{3/2}$', '$4d_{5/2}$']}
dataFiles['BETA'] = {'Data':'Beta', 'File':r'XeF2_Xe4d_beta_all_photon_energies_02102020.dat',
'States':['$\Pi_{1/2} (4d_{3/2})$', '$\Delta_{3/2} (4d_{3/2})$', '$\Sigma_{1/2} (4d_{5/2})$',
'$\Pi_{3/2} (4d_{5/2})$', '$\Delta_{5/2} (4d_{5/2})$']}
dataFiles['BR-All'] = {'Data':'Branching ratios', 'File':r'XeF2_Xe4d_branching_all_photon_energies_all_states_01112020.dat',
'States':['$\Pi_{1/2} (4d_{3/2})$', '$\Delta_{3/2} (4d_{3/2})$', '$\Sigma_{1/2} (4d_{5/2})$',
'$\Pi_{3/2} (4d_{5/2})$', '$\Delta_{5/2} (4d_{5/2})$']}
dataFiles['BR-SO-summed'] = {'Data':'Branching ratios SO summed', 'File':r'XeF2_Xe4d_branching_all_photon_energies_SO_av_01112020.dat',
'States':['$\Pi$', '$\Delta$', '$\Sigma$']}
# Update with ion data
dataFiles2 = {}
dataFiles2['ION-low'] = {'Data':'XS', 'File':r'XeF2_ion_yield_low_energy_cal.txt'}
dataFiles2['ION-high'] = {'Data':'XS', 'File':r'XeF2_ion_yield_high_energy_cal.txt'}
# Update with branching ratios
# dataFilesBR = {}
# dataFilesBR['All'] = {'Data':'Branching ratios', 'File':r'XeF2_Xe4d_branching_all_photon_energies_all_states_01112020.dat',
# 'States':['\pi1/2 (4d3/2)', '\delta3/2 (4d3/2)', '\sigma1/2 (4d5/2)', '\pi3/2 (4d5/2)', '\delta5/2 (4d5/2)']}
# dataFiles['SO-summed'] = {'Data':'Branching ratios SO summed', 'File':r'XeF2_Xe4d_branching_all_photon_energies_SO_av_01112020.dat',
# 'States':['\pi', '\delta', '\sigma']}
# Read data files and convert to Xarray
# 27/10/20 added quick hack to set 2nd array for ion yield data
dataList = []
dataList2 = []
for key in dataFiles:
dataIn = np.loadtxt(dataPathExpt/dataFiles[key]['File'])
# Convert to Xarray
dataXr = xr.DataArray(dataIn[:,1:], dims=('Ehv','State'), coords={'Ehv':dataIn[:,0], 'State':dataFiles[key]['States'][0:dataIn.shape[1]-1]})
dataXr.attrs = dataFiles[key]
dataXr.attrs['dataPath'] = dataPathExpt
dataXr.name = key
dataList.append(dataXr)
# Stack to Xarray
dataExpt = xr.concat(dataList, "XC")
dataExpt['XC'] = list(dataFiles.keys())
for key in dataFiles2:
dataIn = np.loadtxt(dataPathExpt/dataFiles2[key]['File'])
# Convert to Xarray
dataXr = xr.DataArray(dataIn[:,1].squeeze(), dims=('Ehv'), coords={'Ehv':dataIn[:,0]}) # Only 1D datasets in this case
dataXr.attrs = dataFiles2[key]
dataXr.attrs['dataPath'] = dataPathExpt
dataXr.name = key
dataList2.append(dataXr)
dataExpt2 = xr.concat(dataList2, "XC")
dataExpt2['XC'] = list(dataFiles2.keys())
# Quick plot
marker = 'x'
for item in dataExpt:
plt.figure() # Force new figure
item.dropna('State').plot.line(x='Ehv', marker=marker, ls=':') # Include dropna here to remove empty state dims
# Fix axis labels
if item.XC == 'SIGMA':
plt.ylabel('XS/Mb')
elif item.XC == 'BETA':
plt.ylabel(r"$\beta_{LM}$")
else:
plt.ylabel("Branching ratio")
Note dataset labels here - this dataset has XS summed over spin-orbit states (labelled 4d3/2 and 4d5/2), and betas for spin-orbit states.
# Quick plot - Ion yields
marker = 'x'
for item in dataExpt2:
plt.figure() # Force new figure
item.plot.line(x='Ehv') # Include dropna here to remove empty state dims
# Fix axis labels
plt.ylabel("Yield")
Comparison plots...¶
# Compare with computational results
pType = 'SIGMA'
Erange = [50, 250]
pltObj, lTextComp = data.plotGetCroComp(pType=pType, Etype = Etype, Erange = Erange, returnHandles = True)
lText = lTextComp.copy()
# Add expt data - cross-secions
scale = dataExpt.sel({'XC':pType}).max() / data.dataSets['orb21']['XS'].sel({'XC':pType, 'Type':'L'}).max() # Set scaling to match computational data
(dataExpt.sel({'XC':pType})/scale).dropna('State').plot.line(x='Ehv', marker=marker, ls=':');
# Add expt data - ion yields
ionData = 'ION-high'
scale = dataExpt2.sel({'XC':ionData}).max() / data.dataSets['orb21']['XS'].sel({'XC':pType, 'Type':'L'}).max() # Set scaling to match computational data
(dataExpt2.sel({'XC':ionData})/scale).plot.line(x='Ehv');
# Manual legend fix
lText.extend(dataExpt.sel({'XC':pType}).dropna('State').coords['State'].data)
lText.extend([ionData])
plt.legend(lText)
# Fix axis labels
if pType == 'SIGMA':
plt.ylabel('XS/Mb')
plt.title('XS (normalised to computational results)')
else:
plt.ylabel(r"$\beta_{LM}$")
plt.title(r"$\beta_{LM}$")
EDIT 27/10/20: now with ion yield data, which should be taken as the reliable experimental XS result. This doesn't show structure, but is broadly similar to the photoelectron measurements. Structure in photoelectron results to be ignored! Some of the previous comments still apply.
Original notes 13/10/20:
This is interesting... the structures in the computational XS look very similar to the experimental results, but on a much smaller energy scale. I suspect there could be many reasons for this in general (both computational and experimental)... but given that the betas match well (see below), this possibly suggests that this is a signature of multi-electron effects in the bound state.
Some reasons this might be the case:
- The ePS computations are single active electron type calculations, hence do not include any multi-electron/state-mixing effects.
- The good agreement with the betas suggests that the scattering potential, hence continuum wavefunction and angular overlap integrals, are, however, well modelled in the computations.
- This, to me, points to the radial overlap integrals in the computations being quite different to reality, hence anything which strongly influences the radial wavefunction of the bound and/or continuum state. (As a side-note, the bond length used in the calculations could also have an effect here, this was 1.937A, pretty close to the lit value of 1.988A; any other state-mixing, e.g. vibronic couplings, could also play a role of course, although these might generally be expected to affect the angular momentum composition too.)
- Given that the 4d states are mixed, and we also expect multi-electron effects, the bound state seems the most likely candidate for this.
Any suggestions/comments here would be appreciated, since this is rather speculative, and I don't have a good overview of the spectroscopy and general expectations here at all.
04/11/20 - Added branching ratios comparison (rough)...
This shows the computational vs. SO summed experimental branching ratios. Not totally ideal, but should at least show the qualitative trends here. So far it looks roughly OK.
# Compare with computational results
pType = 'BR-SO-summed'
# pltObj, lText = data.plotGetCroComp(pType=pType, Etype = Etype, Erange = Erange, returnHandles = True)
dsXS.swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L').plot.line(x='Ehv');
# Add expt data - cross-secions
# scale = dataExpt.sel({'XC':pType}).max() / data.dataSets['orb21']['XS'].sel({'XC':pType, 'Type':'L'}).max() # Set scaling to match computational data
scale = 1
(dataExpt.sel({'XC':pType})/scale).dropna('State').plot.line(x='Ehv', marker=marker, ls=':');
# Manual legend fix
lText = lTextComp.copy()
lText.extend(dataExpt.sel({'XC':pType}).dropna('State').coords['State'].data)
plt.legend(lText, loc='lower right')
plt.ylabel("Branching ratio")
# Beta comparison plot over orbs
pType = 'BETA'
pltObj, lText = data.plotGetCroComp(pType=pType, Etype = Etype, Erange = Erange, returnHandles = True)
# Add expt data
dataExpt.sel({'XC':pType}).dropna('State').plot.line(x='Ehv', marker=marker, ls=':');
# Manual legend fix
lText.extend(dataExpt.sel({'XC':pType}).dropna('State').coords['State'].data)
plt.legend(lText)
# Fix axis labels
if pType == 'SIGMA':
plt.ylabel('XS/Mb')
plt.title('XS (normalised to computational results)')
else:
plt.ylabel(r"$\beta_{LM}$")
plt.title(r"$\beta_{LM}$")
As previously, applying an energy shift to the data gives a very nice agreement for the betas...
Eshift = 9.0 # NOTE - this breaks energy selection code later, for some reason - float type and/or tolerance in slice?
# Beta comparison plot over orbs
pType = 'BETA'
pltObj, lText = data.plotGetCroComp(pType=pType, Etype = Etype, Erange = Erange, Eshift = Eshift, returnHandles = True)
# Add expt data
dataExpt.sel({'XC':pType}).dropna('State').plot.line(x='Ehv', marker=marker, ls=':');
# Manual legend fix
lText.extend(dataExpt.sel({'XC':pType}).dropna('State').coords['State'].data)
plt.legend(lText)
# Fix axis labels
if pType == 'SIGMA':
plt.ylabel('XS/Mb')
else:
plt.ylabel(r"$\beta_{LM}$")
Spin-orbit calculations¶
Formalism¶
For ion SO case only, should be very similar to "molecular" form from old (old) work, so use this as a paradigm. In terms of ion electronic state only, using Hund's case b/c notation (essentially identical if one neglects rotational ang. mom.):
\begin{equation} C^{SO}(L,J,S)=\left(\begin{array}{ccc} L & J & S\\ M_{L} & M_{J} & M_{S} \end{array}\right)\left(\begin{array}{ccc} L & J & S\\ \Lambda & \Omega & \Sigma \end{array}\right) \end{equation}And coherent sum:
\begin{equation} \Xi(L,J,S)=\sum_{all\,projections}C^{SO}(L,J,S)C^{SO}(L',J',S')=\sum_{all\,projections}\left(\begin{array}{ccc} L & J & S\\ M_{L} & M_{J} & M_{S} \end{array}\right)\left(\begin{array}{ccc} L & J & S\\ \Lambda & \Omega & \Sigma \end{array}\right)\left(\begin{array}{ccc} L' & J' & S'\\ M'_{L} & M'_{J} & M'_{S} \end{array}\right)\left(\begin{array}{ccc} L' & J' & S'\\ \Lambda' & \Omega' & \Sigma' \end{array}\right) \end{equation}All states will be modulated by coherent sum - if this is summed over all projection terms then it collapses to a $6j$ (need to check phases carefully here however!). If a single set of $(L,J,S)$ are assumed (i.e. states are resolved), $L=L',\,J=J',\,S=S'$. Assuming all (lab frame) $M$ are equally populated, only the MF term will affect things, so we'll only need the square of the term with $(\Lambda,\Omega,\Sigma)$ in the present case:
\begin{eqnarray} \Xi^{MF}(L,J,S) & = & \sum_{all\,unresolved}\left(\begin{array}{ccc} L & J & S\\ \Lambda & \Omega & \Sigma \end{array}\right)\left(\begin{array}{ccc} L & J & S\\ \Lambda' & \Omega' & \Sigma' \end{array}\right)\\ & = & \left(\begin{array}{ccc} L & J & S\\ \Lambda & \Omega & \Sigma \end{array}\right)^{2} \end{eqnarray}Where the first line is applicable in the case of unresolved states, and the second if all QNs are defined.
TODO: derive this properly... likely missing some sign conventions/phases here, although will fall out in the numerics if all +- term combinations are included.
Refs (see also refs. therein):
Hockett, Paul, and Katharine L Reid. “Complete Determination of the Photoionization Dynamics of a Polyatomic Molecule. II. Determination of Radial Dipole Matrix Elements and Phases from Experimental Photoelectron Angular Distributions from A1Au Acetylene.” The Journal of Chemical Physics 127, no. 15 (October 2007): 154308. https://doi.org/10.1063/1.2790443.
Hockett, Paul. “Photoionization Dynamics of Polyatomic Molecules.” PhD Thesis, University of Nottingham, 2009. http://eprints.nottingham.ac.uk/10857/.
Application to $XeF_{2}(4d^{-1})$¶
In this case we set the following for the spin-orbit splitting in the ion (TBC!):
- $L, S$ corresponds to the $4d$ unpaired electon,
- $J = L\pm S$
- Hence the overall term is $^{2s+1}L_{J}=^{2}D_{3/2,5/2}$
- $\Lambda$ corresponds to the ligand-field split components (equivalently, the components of the 4d orbital), $\Lambda = 0, 1, 2$, corresponding to $\Sigma, \Pi, \Delta$ components or, equivalently, the A1G/SG, E1G/PG, E2G/DG ab initio states respectively.
Tests¶
from epsproc.geomFunc.geomUtils import genllLList
from epsproc.geomFunc.geomCalc import w3jTable
# Gen QNs for specific (L,S) case
L = 2
S = 0.5
Llist = np.array([[L,L+S,S], [L, L, S], [L,L-S,S]]) # Note this needs to be 2D array in current form of function
QNs = genllLList(Llist, uniqueFlag = False)
# Then calc 3js....
backend = 'sympy'
form = 'xdaLM' # '2d' # 'xdaLM' # xds
nonzeroFlag = True
dlist = ['L', 'J', 'S', 'Lambda', 'Omega', 'Sigma'] # Set dims for reference
thrj = w3jTable(QNs = QNs, nonzeroFlag = nonzeroFlag, form = form, dlist = dlist, backend = backend)
# And primed terms (will be identical at this point, but set dims for multiplication later)
thrjP = w3jTable(QNs = QNs, nonzeroFlag = nonzeroFlag, form = form, dlist = [item + 'p' for item in dlist], backend = backend)
# BASIC case of single 3j term
pdTable, _ = ep.multiDimXrToPD(thrj, colDims = 'Lambda')
pdTable
This seems OK... but has more allowed states than the experimental assignments, in particular $J=3/2,\Sigma_1/2$ is non-zero, and there are additional terms in $\Omega$. This all seems to be correct for the 3j term alone, so I'm presumably just missing some other effects here (e.g. parity restriction), and the $\Omega$ terms maybe degenerate (?). It could also be an issue with the choice of Hund's case coupling here... I need to look into the spectroscopy a bit more here to confirm.
For now, we'll assume the best and see how the results look...
# Reformate by comsolidating +/- terms as modulation for XS and square
thrjUS = thrj.unstack('mSet').fillna(0)
# thrjXSmod = 2*(thrjUS.sum('Sigma').where((thrjUS['Lambda']>-0.5)&(thrjUS['Omega']<0.5), drop=True))**2 # Note Lambda & Omega anti-phase!
thrjXSmod = 2*(thrjUS.sum('Sigma').where((thrjUS['Lambda']>-0.5)&(thrjUS['Omega']<0.5), drop=True))**2 # Note Lambda & Omega anti-phase!
thrjXSmod['Omega'] = np.abs(thrjXSmod['Omega'])
pdTable, _ = ep.multiDimXrToPD(thrjXSmod, colDims = 'Lambda')
pdTable
For the coherent case, allow all terms in projection pairs (unprimed & primed), then sum and/or subselect as required.
# Terms for coherent case
thrjSumCoherent = (thrj).unstack('mSet').fillna(0) * (thrjP).unstack('mSet').fillna(0)
# Sum over spin only, specify Omega = Omega', and Lambda = Lambda'
thrjSumCoherent = thrjSumCoherent.sum(['Sigmap', 'Sigma']).where(thrjSumCoherent['Omega']==thrjSumCoherent['Omegap']).where(thrjSumCoherent['Lambda']==thrjSumCoherent['Lambdap'])
# Drop (sum) redundant coords & reindex
thrjXSmodCoherent = 2*thrjSumCoherent.sum(['Omegap','Lambdap']).where((thrjUS['Lambda']>-0.5)&(thrjUS['Omega']<0.5), drop=True)
thrjXSmodCoherent['Omega'] = np.abs(thrjXSmodCoherent['Omega'])
pdTable, _ = ep.multiDimXrToPD(thrjXSmodCoherent, colDims = 'Lambda') #
pdTable
Same result here, since the sum over $\Sigma,\Sigma'$ is just a degeneracy factor, and other terms are set to be equal (no additional coherences), which shows the numerics are OK!
If we assume that $\Omega$ remains coherent (unresolved)...
# Terms for coherent case
thrjSumCoherent = (thrj).unstack('mSet').fillna(0) * (thrjP).unstack('mSet').fillna(0)
# Sum over spin only, specify Lambda = Lambda'
thrjSumCoherent = thrjSumCoherent.sum(['Sigmap', 'Sigma']).where(thrjSumCoherent['Lambda']==thrjSumCoherent['Lambdap'])
# pdTable, _ = ep.multiDimXrToPD(thrjSumCoherent, colDims = 'Lambda') #
# Drop (sum) redundant coords & reindex
thrjXSmodCoherent = 2*thrjSumCoherent.sum(['Omegap','Lambdap']).where((thrjUS['Lambda']>-0.5)&(thrjUS['Omega']<0.5), drop=True)
thrjXSmodCoherent['Omega'] = np.abs(thrjXSmodCoherent['Omega'])
pdTable, _ = ep.multiDimXrToPD(thrjXSmodCoherent, colDims = 'Lambda') #
pdTable
Here the $\Sigma_{1/2}$ term goes to zero for $J=5/2$... this might be genuine, or simply a phase issue in the summation! There are also negative values present, which likely points to the same issue.
If we treat this as per Hund's case (c), and assume it is $\Lambda$ which is unresolved/undefined,
# Terms for coherent case
thrjSumCoherent = (thrj).unstack('mSet').fillna(0) * (thrjP).unstack('mSet').fillna(0)
# Sum over spin only, specify Omega = Omega'
thrjSumCoherent = thrjSumCoherent.sum(['Sigmap', 'Sigma']).where(thrjSumCoherent['Omega']==thrjSumCoherent['Omegap'])
# pdTable, _ = ep.multiDimXrToPD(thrjSumCoherent, colDims = 'Lambda') #
# Drop (sum) redundant coords & reindex
thrjXSmodCoherent = 2*thrjSumCoherent.sum(['Omegap','Lambdap']).where((thrjUS['Lambda']>-0.5)&(thrjUS['Omega']<0.5), drop=True)
thrjXSmodCoherent['Omega'] = np.abs(thrjXSmodCoherent['Omega'])
pdTable, _ = ep.multiDimXrToPD(thrjXSmodCoherent, colDims = 'Lambda') #
pdTable
SO Branching ratios¶
Assuming that the above is "correct" the branching ratios should follow fairly directly...
First, define a basic function to multiply by 3j term(s) and set branching ratios. Then test for some different cases...
# Test SO multiplication terms
# 11/11/20 - converted to function
def simulateSOBR(dsXS, thrjSOTerm, BRtype = 'All', coupling = 'Lambda', stateLabels = ['$\Sigma$', '$\Pi$', '$\Delta$']):
# Set Lambda terms
# dsXS['Lambda'] = dsXS['channel']
# dsXS['Lambda'] = [0.0, 1.0, 2.0]
dsXSO = dsXS.assign_coords({coupling:('channel',thrjSOTerm[coupling])}).swap_dims({'channel':coupling}) # Assign Lambda for multiplication
dsXSO = dsXSO.assign_coords({f'{coupling}-Labels':(coupling, np.asarray(stateLabels)[dsXSO[coupling].astype(int).data])}) # Assign labels, based on int values as indexers (CRUDE!)
dsXSO = dsXSO * thrjSOTerm
# dsXS
# Convert to branching ratios (total)
if BRtype == 'All':
# dsXSO['total'] = dsXSO.sum(['Lambda', 'Omega', 'lSet']) # Sum over channels
dsXSO['total'] = dsXSO.sum(thrjSOTerm.dims)
elif BRtype == 'J':
dsXSO['total'] = dsXSO.sum(['Lambda', 'Omega']) # Sum over channels, except J - GIVES same result?
elif BRtype == 'None':
dsXSO['total'] = 1.0
else:
print(f'BRtype={BRtype} not recognised.')
return
# # Set branching ratios
dsXSO = dsXSO/dsXSO['total']
return dsXSO
Assuming $\Lambda$ good QN¶
Label levels by $\Lambda$, and use all $\Omega$ components in summation.
# For the state-resolved case
dsXSO = simulateSOBR(dsXS, thrjXSmod, BRtype = 'J')
# Plot
Etype = 'Ehv'
Erange = [75, 250]
dsXSO.swap_dims({'Eke':'Ehv', 'Lambda':'Lambda-Labels'}).sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L').unstack().squeeze().plot.line(x='Ehv', col = 'Omega', row = 'J');
# dsXSO.swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L').unstack().squeeze().plot.line(x='Ehv', col = 'Omega', row = 'J');
# Summed over Omega
# Note J states look essentially identical here, as they should in the absence of any other modulating factors!
dsXSO.swap_dims({'Eke':'Ehv', 'Lambda':'Lambda-Labels'}).sum('Omega').sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L').unstack().squeeze().plot.line(x='Ehv', row = 'J');
# Compare with computational results
pType = 'BR-All'
# pltObj, lText = data.plotGetCroComp(pType=pType, Etype = Etype, Erange = Erange, returnHandles = True)
# Loop over allowed Omega terms and add to plot
# for omega in dsXSO.Omega:
# J=2.5
# (2*dsXSO).swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L', J=J, Omega=omega).squeeze().plot.line(x='Ehv');
# Loop over J, sum Omega
# for J in dsXSO.J:
# (dsXSO).sum('Lambda').swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L', J=J.item()).squeeze().plot.line(x='Ehv');
# Selected term(s) only
J=1.5
dsPlot = (dsXSO).swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L', J=J).sum('Omega').squeeze()
dsPlot.plot.line(x='Ehv');
# Add expt data - cross-secions
# scale = dataExpt.sel({'XC':pType}).max() / data.dataSets['orb21']['XS'].sel({'XC':pType, 'Type':'L'}).max() # Set scaling to match computational data
scale = 1
# All data
# (dataExpt.sel({'XC':pType})/scale).dropna('State').plot.line(x='Ehv', marker=marker, ls=':');
# Selected states only
statesPlot = ['$\\Delta_{3/2} (4d_{3/2})$', '$\\Pi_{1/2} (4d_{3/2})$']
(dataExpt.sel({'XC':pType, 'State':statesPlot})/scale).plot.line(x='Ehv', marker=marker, ls=':');
# Manual legend fix
lText = list(dsPlot['Lambda-Labels'].data)
lText.extend(statesPlot)
# lText.extend(dataExpt.sel({'XC':pType}).dropna('State').coords['State'].data)
plt.legend(lText, loc='lower right')
plt.ylabel("Branching ratio")
plt.title(f'Simulated branching ratios, J={J}, Hund\'s case (c)')
# Compare with computational results
pType = 'BR-All'
# pltObj, lText = data.plotGetCroComp(pType=pType, Etype = Etype, Erange = Erange, returnHandles = True)
# Loop over allowed Omega terms and add to plot
# for omega in dsXSO.Omega:
# J=2.5
# (2*dsXSO).swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L', J=J, Omega=omega).squeeze().plot.line(x='Ehv');
# Loop over J, sum Omega
# for J in dsXSO.J:
# (dsXSO).sum('Lambda').swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L', J=J.item()).squeeze().plot.line(x='Ehv');
# Selected term(s) only
J=2.5
(dsXSO).swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L', J=J).sum('Omega').squeeze().plot.line(x='Ehv');
# Add expt data - cross-secions
# scale = dataExpt.sel({'XC':pType}).max() / data.dataSets['orb21']['XS'].sel({'XC':pType, 'Type':'L'}).max() # Set scaling to match computational data
scale = 1
# All data
# (dataExpt.sel({'XC':pType})/scale).dropna('State').plot.line(x='Ehv', marker=marker, ls=':');
# Selected states only
statesPlot = ['$\\Delta_{5/2} (4d_{5/2})$', '$\\Pi_{3/2} (4d_{5/2})$', '$\\Sigma_{1/2} (4d_{5/2})$']
(dataExpt.sel({'XC':pType, 'State':statesPlot})/scale).plot.line(x='Ehv', marker=marker, ls=':');
# Manual legend fix
lText = list(dsPlot['Lambda-Labels'].data)
lText.extend(statesPlot)
# lText.extend(dataExpt.sel({'XC':pType}).dropna('State').coords['State'].data)
plt.legend(lText, loc='lower right')
plt.ylabel("Branching ratio")
plt.title(f'Simulated branching ratios, J={J}, Hund\'s case (c)')
This does not look very convincing... too many states, and not a good match to the experimental BRs.
Assuming $\Omega$ good QN (should be correct for Hund's case c), and that it labels states...¶
Label levels by $\Omega$, and use all $\Lambda$ components in summation.
# For the state-resolved case
# *** WITH Omega & Lambda switched (Hund's case c)
# This seems to give better results, but may still be incorrect - should just subselect terms instead of sum?
dsXSO = simulateSOBR(dsXS, thrjXSmod, BRtype = 'J', coupling = 'Omega')
# Plot
Etype = 'Ehv'
Erange = [75, 250]
dsXSO.swap_dims({'Eke':'Ehv', 'Omega':'Omega-Labels'}).sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L').unstack().squeeze().plot.line(x='Ehv', col = 'Lambda', row = 'J');
# dsXSO.swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L').unstack().squeeze().plot.line(x='Ehv', col = 'Lambda', row = 'J');
# Summed over Lambda
# Note J states look essentially identical here, as they should in the absence of any other modulating factors!
dsXSO.swap_dims({'Eke':'Ehv', 'Omega':'Omega-Labels'}).sum('Lambda').sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L').unstack().squeeze().plot.line(x='Ehv', row = 'J');
# Compare with computational results
pType = 'BR-All'
# pltObj, lText = data.plotGetCroComp(pType=pType, Etype = Etype, Erange = Erange, returnHandles = True)
# Loop over allowed Omega terms and add to plot
# for omega in dsXSO.Omega:
# J=2.5
# (2*dsXSO).swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L', J=J, Omega=omega).squeeze().plot.line(x='Ehv');
# Loop over J, sum Omega
# for J in dsXSO.J:
# (dsXSO).sum('Lambda').swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L', J=J.item()).squeeze().plot.line(x='Ehv');
# Selected term(s) only
J=1.5
dsPlot = (dsXSO).swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L', J=J).sum('Lambda').squeeze()
dsPlot.plot.line(x='Ehv');
# Add expt data - cross-secions
# scale = dataExpt.sel({'XC':pType}).max() / data.dataSets['orb21']['XS'].sel({'XC':pType, 'Type':'L'}).max() # Set scaling to match computational data
scale = 1
# All data
# (dataExpt.sel({'XC':pType})/scale).dropna('State').plot.line(x='Ehv', marker=marker, ls=':');
# Selected states only
statesPlot = ['$\\Delta_{3/2} (4d_{3/2})$', '$\\Pi_{1/2} (4d_{3/2})$']
(dataExpt.sel({'XC':pType, 'State':statesPlot})/scale).plot.line(x='Ehv', marker=marker, ls=':');
# Manual legend fix
lText = list(dsPlot['Omega-Labels'].data)
lText.extend(statesPlot)
# lText.extend(dataExpt.sel({'XC':pType}).dropna('State').coords['State'].data)
plt.legend(lText, loc='lower right')
plt.ylabel("Branching ratio")
plt.title(f'Simulated branching ratios, J={J}, Hund\'s case (c)')
# Compare with computational results
pType = 'BR-All'
# pltObj, lText = data.plotGetCroComp(pType=pType, Etype = Etype, Erange = Erange, returnHandles = True)
# Loop over allowed Omega terms and add to plot
# for omega in dsXSO.Omega:
# J=2.5
# (2*dsXSO).swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L', J=J, Omega=omega).squeeze().plot.line(x='Ehv');
# Loop over J, sum Omega
# for J in dsXSO.J:
# (dsXSO).sum('Lambda').swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L', J=J.item()).squeeze().plot.line(x='Ehv');
# Selected term(s) only
J=2.5
(dsXSO).swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L', J=J).sum('Lambda').squeeze().plot.line(x='Ehv');
# Add expt data - cross-secions
# scale = dataExpt.sel({'XC':pType}).max() / data.dataSets['orb21']['XS'].sel({'XC':pType, 'Type':'L'}).max() # Set scaling to match computational data
scale = 1
# All data
# (dataExpt.sel({'XC':pType})/scale).dropna('State').plot.line(x='Ehv', marker=marker, ls=':');
# Selected states only
statesPlot = ['$\\Delta_{5/2} (4d_{5/2})$', '$\\Pi_{3/2} (4d_{5/2})$', '$\\Sigma_{1/2} (4d_{5/2})$']
(dataExpt.sel({'XC':pType, 'State':statesPlot})/scale).plot.line(x='Ehv', marker=marker, ls=':');
# Manual legend fix
lText = list(dsPlot['Omega-Labels'].data)
lText.extend(statesPlot)
# lText.extend(dataExpt.sel({'XC':pType}).dropna('State').coords['State'].data)
plt.legend(lText, loc='lower right')
plt.ylabel("Branching ratio")
plt.title(f'Simulated branching ratios, J={J}, Hund\'s case (c)')
This looks quite reasonable, with the same number (although different ordering) of states assigned experimentally. Could still be coupling case problem? Also looks to be a relative phase issue/flip in J=3/2 case?
State selected Assuming both $\Lambda$ and $\Omega$ are defined... as per experimental assignments (?)¶
Label levels by ($\Omega$, $\Lambda$), with no additional summation over terms. Q: is this justified/realistic?
# Selected term(s) only - LOOKS GOOD FOR coupling='Omega' case, but opposite phases for coupling='Lambda'
dsXSO = simulateSOBR(dsXS, thrjXSmod, BRtype = 'None', coupling='Omega') # This looks good
# dsXSO = simulateSOBR(dsXS, np.abs(thrjXSmodCoherent), BRtype = 'None', coupling='Lambda') # Testing coherent case
# Unnorm
# dsXSO *= dsXSO['total']
# Selected term(s) only
J=1.5
OLset = [[2,1.5], [1,0.5]] # [Omega, Lambda] pairs to match expt.
# J=2.5
# OLset = [[2,2.5], [1,1.5], [0,0.5]] # [Omega, Lambda] pairs to match expt.
# subset = xr.DataArray()
subset = []
for OL in OLset:
subset.append((dsXSO).swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L', J=J, Lambda=OL[0], Omega=OL[1]))
subsetXS = xr.concat(subset, dim='Omega').squeeze()
# Renorm
subsetXS = subsetXS/subsetXS.sum('Omega')
subsetXS.plot.line(x=Etype)
# Selected states only
if J == 1.5:
statesPlot = ['$\\Delta_{3/2} (4d_{3/2})$', '$\\Pi_{1/2} (4d_{3/2})$']
else:
statesPlot = ['$\\Delta_{5/2} (4d_{5/2})$', '$\\Pi_{3/2} (4d_{5/2})$', '$\\Sigma_{1/2} (4d_{5/2})$']
(dataExpt.sel({'XC':pType, 'State':statesPlot})/scale).plot.line(x='Ehv', marker=marker, ls=':');
# Manual legend fix
# lText = list(subsetXS['Omega'].data)
lText = list(OLset)
# lText = list(subsetXS['Lambda'].data)
lText.extend(statesPlot)
# lText.extend(dataExpt.sel({'XC':pType}).dropna('State').coords['State'].data)
plt.legend(lText, loc='lower right')
plt.ylabel("Branching ratio")
plt.title(f'Simulated branching ratios, J={J}, Hund\'s case (c), with subselection (no summation)')
# Selected term(s) only
# Selected term(s) only - LOOKS GOOD FOR coupling='Omega' case, but opposite phases for coupling='Lambda'
# For J = 5/2 'Lambda' case is possibly better, at least for state ordering? DO NOT UNDERSTAND THIS YET!
dsXSO = simulateSOBR(dsXS, thrjXSmod, BRtype = 'None', coupling='Omega')
# Unnorm
# dsXSO *= dsXSO['total']
# Selected term(s) only
# J=1.5
# OLset = [[2,1.5], [1,0.5]] # [Omega, Lambda] pairs to match expt.
J=2.5
OLset = [[2,2.5], [1,1.5], [0,0.5]] # [Omega, Lambda] pairs to match expt.
# subset = xr.DataArray()
subset = []
for OL in OLset:
subset.append((dsXSO).swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L', J=J, Lambda=OL[0], Omega=OL[1]))
subsetXS = xr.concat(subset, dim='Omega').squeeze()
# Renorm
subsetXS = subsetXS/subsetXS.sum('Omega')
subsetXS.plot.line(x=Etype)
# Selected states only
if J == 1.5:
statesPlot = ['$\\Delta_{3/2} (4d_{3/2})$', '$\\Pi_{1/2} (4d_{3/2})$']
else:
statesPlot = ['$\\Delta_{5/2} (4d_{5/2})$', '$\\Pi_{3/2} (4d_{5/2})$', '$\\Sigma_{1/2} (4d_{5/2})$']
(dataExpt.sel({'XC':pType, 'State':statesPlot})/scale).plot.line(x='Ehv', marker=marker, ls=':');
# Manual legend fix
# lText = list(subsetXS['Omega'].data)
lText = list(OLset)
# lText = list(subsetXS['Lambda'].data)
lText.extend(statesPlot)
# lText.extend(dataExpt.sel({'XC':pType}).dropna('State').coords['State'].data)
plt.legend(lText, loc='lower right')
plt.ylabel("Branching ratio")
plt.title(f'Simulated branching ratios, J={J}, Hund\'s case (c), with subselection (no summation)')
Not too bad for Hund's case (c), plus subselecting to match only the states assigned experimentally. This looks really nice for J=3/2, and not too bad for J=5/2, although I get different state ordering in that case (treating as Hund's case (b) possibly gives a better result in this case, but is that justified?).
q's
- Is the formalism correct?
- Are these numerics correct... may be issue with final subselected results due to choice of leading dim in multiplication?
- SPECTROSCOPY?
- Experimental assignments OK?
- Experimental data uncertainties, might expect more bleed between partially-resolved J=5/2 features.
Addendum¶
Molecular orbitals¶
SEE: https://epsproc.readthedocs.io/en/dev/methods/ePSproc_orbPlot_tests_130520.html
from epsproc.vol import orbPlot
# chemPath = r'/home/femtolab/python/chem/tools/chemlab/chemlab' # Linux dev machine
chemPath = r'D:\temp\chemlab\chemlab' # Win dev machine
out = orbPlot.importChemLabQC(chemPath = chemPath)
# Test class
# filename = r"/home/femtolab/ePS/XeF2/electronic_structure/xef2_SPKrATZP_rel_geom.log" # XeF2 Gamess file OK
filename = r"D:\projects\ePolyScat\xef2\electronic_structure\xef2_SPKrATZP_rel_geom.log"
mo = orbPlot.molOrbPlotter(chemPath = chemPath, fileIn = filename)
# Set orbs to plot - currently handles list of strs only!
orbs = [21,22,23,24,25]
orbList = [str(item) for item in orbs]
[mo.calcOrb(orbN = orbN) for orbN in orbs];
# Plot orbs as independent figures
[mo.plotOrb(orbN = [item], interactive = False, showAtoms = False) for item in orbList];
This shows orbitals 21 - 25, the full set of 4d components. Note that (22,23) and (24,25) are degenerate pairs, and are treated as one orbital (with occ = 4) in ePolyScat, as noted previously.
Gauge results comparison¶
This shows a big change in absolute X-section between lenght (L) and velocity (V) gauge calculations (~order of magnitue), but the feature are appoximately the same in both cases (apologies for the poor plots here!).
data.plotGetCro(pType='SIGMA', Etype=Etype, Erange=Erange)
# Branching ratios
dsXS.swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}).plot.line(x='Ehv', col='Type');
Beta parameters are pretty much idendical for both gauges, with a broad, shape-resonance type feature at lower photon energies, and lots of small-scale resonance structure to higher photon energies.
data.plotGetCro(pType='BETA', Etype=Etype, Erange=Erange)
XS and $\beta$ per orbital & symmetry¶
Interactive^ plots per orbital & symmetry. In these plots the solid lines show the 'mixed' gauge results, and dashed lines show length and velocity gauges as bounds on the values.
^ Use controls in top right of plot. Also data sets can be turned on/off using the legend entries.
data.plotGetCro(pType='SIGMA', Etype = Etype, Erange = Erange, backend = 'hv')
Matrix elements¶
data.lmPlot(dataType = 'matE', Erange = Erange, Etype = Etype, thres = 0.1, logFlag = False, selDims = {'Type':'L'}, fillna = True)
Versions¶
import scooby
scooby.Report(additional=['epsproc', 'xarray', 'jupyter'])